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Analysis of an asymptotic preserving scheme for linear kinetic equations in the 

diffusion limit 



Jian-Guo Liv0, Luc Mieussen^ 

O - 

O ■ Abstract. We present a mathematical analysis of the asymptotic preserving scheme 

proposed in [M. Lemou and L. Mieussens, SI AM J. Sci. Comput., 31, pp. 334-368, 2008] 

I for linear transport equations in kinetic and diffusive regimes. We prove that the scheme 

O [ is uniformly stable and accurate with respect to the mean free path of the particles. This 

■ property is satisfied under an explicitly given CFL condition. This condition tends to a 

parabolic CFL condition for small mean free paths, and is close to a convection CFL condition 

^ ! for large mean free paths. Our analysis is based on very simple energy estimates. 
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oo : 1 Introduction 

^ ■ Particle systems are often described at the microscopic level by kinetic models (neutron 
^ . transport, radiative transfer, electrons in semi-conductors, or rarefied gas dynamics). Sim- 
ulating such systems by using kinetic models can be computationally very expensive, but 
^ I modern super computers now enable realistic simulations. When the mean free path of the 
particles is very small as compared to the (macroscopic) size of the computational domain, 
the kinetic model can be very well approximated by a much simpler macroscopic model (dif- 
fusion equation, Rosseland approximation, Euler and Navier-Stokes equations), that can be 
5^ i numerically solved much faster. 

However, there are many cases where the ratio mean free path/macroscopic size (the 
so-called "Knudsen number" in rarefied gas dynamics, denoted by e in this paper) is not 
constant: depending on the geometry of the boundaries, or on the boundary conditions, this 
ratio may vary in time, and in space. In such multiscale situations, usual kinetic solvers are 
often useless: for stability and accuracy reasons, they must resolve the microscopic scales, 
which is computationally to expensive in "fluid" zones (where e is small). By contrast, 
macroscopic solvers are faster but may be inaccurate in "kinetic" zones (where e is large). 
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This is why many people have been working for more than 20 years on a kind of multiscale 
kinetic schemes: the asymptotic-preserving (AP) schemes. Such schemes are uniformly 
stable with respect to e (thus their computational complexity does not depend on e), and 
are consistent with the macroscopic model when e goes to (the limit of the scheme is a 
scheme for the macroscopic model). 

Up to our knowledge, AP schemes have first been studied (for steady problems) in neutron 
transport by Larsen, Morel and Miller p2], Larsen and Morel [21], and then by Jin and 
Levermore [lOl |TT]. For unstationary problems, the difficulty is the time stiffness due to 
the collision operator. To avoid the use of expensive fully implicit schemes, two classes of 
semi-implicit time discretizations have been proposed by Klar |16] and Jin, Pareschi and 
Toscani [15] (see preliminary works in [TH [9] and extensions in [T3l [T2| [25l [TTf [18] . Similar 
ideas have also been used by Buet et al. in [1]. 

In [23] , Lemou and Mieussens have proposed a new AP scheme based on the micro- macro 
decomposition of the distribution function into microscopic and macroscopic components 
(similar schemes have also been proposed by Klar and Schmeiser [T^], and more recently 
by Carrillo et al. [3 [6]). A coupled system of equations is obtained for these two compo- 
nents without any linearity assumption. The decomposition only uses basic properties of 
the collision operator that are common to most of kinetic equations (namely conservation 
and equilibrium properties). Then this system is solved with a suitable time semi- implicit 
discretization and space finite differences on staggered grids. While almost all the schemes 
mentioned before are based on very similar ideas, the approach proposed in [23] has been 
shown to be very general, since it applies to kinetic equations for both diffusion and hy- 
drodynamic regimes (for the diffusion regime, see [23j for linear transport equations and [3j 
for the non-linear Kac equation, for the hydrodynamic regime, see [2] for the Boltzmann 
equation). We also mention the work of Degond, Liu and Mieussens [7] who proposed a 
similar approach (micro-macro decomposition) to design macroscopic diffusion models with 
kinetic upscalings: this approach also leads to AP schemes, at least for a semi-discrete time 
discretization (no AP space discretization was studied in this paper). 

While many different schemes have been proposed in the past few years, it appears that 
the rigorous proof of their AP property looks rather difficult and is seldom investigated. Up 
to our knowledge, there are only two papers on this subject. Klar and Unterreiter [20] have 
proved that a scheme similar to that of [16] and [15] for the linear transport equation is 
uniformly stable. However, their proof is based on a von Neumann analysis, and hence is re- 
stricted to a one dimensional equation with constant coefficients in a periodic domain. Gosse 
and Toscani have proposed in [S] an AP scheme based on a rather different idea (discretiza- 
tion of the collision term as a non-conservative product and use of well-balanced Godunov 
schemes): they have been able to prove a uniform stability property, still in the linear case, 
but also a strong positivity property. However, their scheme is based on techniques (like 
approximation of the steady solution) that are difficult to generalize to other equations. 

In this paper, we propose a very simple stability proof for the AP scheme of [23J and we 
exhibit an explicit CFL condition. This condition is uniform with respect to e, and gives 
a standard parabolic CFL condition At = 0(Aa:^) when e goes to 0. While a stability 
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property was already proved in [23] for this scheme, this was only for a simple two-velocity 
model (telegraph equation), by using von Neumann analysis. Here, our proof is based on 
energy estimates that are more general than von Neumann analysis, and hence is valid for 
one-dimensional linear equations with non-constant coefficients, continuous velocity variable, 
in the whole space. By the same technique, we are also able to prove uniform error estimates. 

Our paper is organized as follows. In section [21 we introduce a general linear equation, 
we present its discretization by the AP scheme, and we give the main features of this scheme. 
Then, in section [3l we give our main stability result and its proof. The error estimates are 
given in section [H 

2 An AP scheme for the linear transport equation 

Linear transport equation is a model for the evolution of particles in some medium (neutron 
transport, linear radiative transfer, . . . ). Generally, this model reads, in scaled variables 



where (p{t, r, Q) is the number density of particles in the position-direction phase space that 
depends on time t, position r = {x,y,z) G M^, and angular direction of propagation of 
particles Q & S"^. Moreover, a is the total cross section, a a is the absorption cross section, 
and S is an internal source of particles, which is independent of fl. The linear operator L 
models the scattering of the particles by the medium and acts only on the angular dependence 
of (I). This simple model does not allow for particles of possibly different energy (or frequency); 
it is called "one-group" or "monoenergetic" equation. The parameter e is a scale factor that 
measures the ratio between a typical microscopic length (the mean free path of a particle, 
for instance) to a typical macroscopic length (the size of the computational domain, for 
instance). See [22] for details. 

In this paper, we consider this one-group equation in the slab geometry: we assume that 
depends only on the slab axis variable x G M. Then it can be shown that the average 
of (j) with respect to the {y,z) cosine directions of Q, denoted by f{t,x,v), satisfies the 
one-dimensional equation 



where v G [—1, 1] is the x cosine direction of fl. Att = 0, we have the initial data /(O, x, v) = 
/°(a;,f). We assume that the cross sections satisfy the inequalities < 0"^ < cr^x) < (Tm 
and < (Ta{x) < a am for every x. We do not consider boundary conditions in this paper, 
and hence a; G M. 

The linear operator L is given by 



edtcf) + ^ ■ Vr<p = —Lcf) — eaAcj) + eS, 
e 



edtf + vd^f = -Lf - eaAf + eS, 





(2) 
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where the kernel s is such that < Sm < s{v,v') < sm for every v,v' in [—1, 1]. We also 
assume that s satisfies J^-^^ s{v , v') dv' = 1, and that it is symmetric: s{v,v') = s{v,v'). For 

the sequel, it is useful to define the operator [.] such that [0] = | J^^ (j){y) dv is the average 
of every velocity dependent function 0. With these assumptions, it is standard [T] to state 
the following properties of L: 

Proposition 2.1. • [Lcf)] = for every (p in L^([— 1, 1]) 

• the null space of L is Af{L) = {(/)= [(f)]} (constant functions) 

• the rank of L is n{L) = U^{L) = {(j) s.t [0] = 0} 

• L is non-positive self-adjoint in 1, 1]) and we have 

[(t>L(t)] < -2sm [(^'] (3) 

for every (p G Af^{L) 

• L admits a pseudo-inverse from //^{L) onto N'^{L), denoted by 

• the orthogonal projection from L'^{[— 1,1]) onto J\f-^{L) is [.] 

When e becomes small (the "diffusion" regime), it is well known that the solution / of ([I]) 
tends to its own average density p = [f], which is a solution of the asymptotic diffusion limit 

dtp - d^Kd^p = -Gap + (4) 

where the diffusion coefficient is k{x) = An asymptotic preserving scheme for the 

linear kinetic equation ([1]) is a numerical scheme that discretizes ([T]) in such a way that it 
leads to a correct discretization of the diffusion limit (jll) when e is small. 

Now, we summarize the results obtained in [23j. By using the micro-macro decomposition 
f = p -\- eg, where p = [/] and g is such that [^r] = 0, we derived the micro-macro model 
for ([1]) that reads 

dtp + [vg] = -Gap + S, (5a) 
dtg + -(/ - []){vd^g) = -^Lg - \vd^p, (5b) 

with initial data p° = [f^] and eg^ = f^ — P^- This system if formally equivalent to ([1]). 

Then, we proposed the following numerical scheme for this system. We choose a time 
step At and times in = nAt, and two staggered grids of step Ax and nodes Xi = iAx and 
= {i-\-^)Ax. We use the approximated values ^ p(tn, Xi) and g^,i{v) ~ g{tn, , v). 
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and the scheme reads 



At 



+ 



* 2 



Ax 



At 



(6a) 
(6b) 



1 p-ri _ n 
^ Pi+1 Pi 



where = Note that as in the continuous case, this scheme preserves the zero average 

of g: 



Proposition 2.2. If the initial data satisfies 
i, we have 



0. 



for every i, then for every n and 



Proof. Apply the average operator [.] to ([6b]): since we have [/ 
(proposition 12. ip . and [v] = (obvious), this relation yields 
the result. 



(7) 

(obvious), [L] = 
0, which gives 
□ 



n 



In scheme ([6]), the upwind discretization of (/ — [.]){vdxg) is to insure stability in the 
kinetic regime, while the centered approximations of [vg] and vd^p are to capture the 



diffusion limit. Indeed, it is clear that when e goes to 0, we have from (l6bl) g^^l 



Ax 



+ 0(e). Consequently, the flux of g^^l is 



1 1 Pi+1 Pi 



Ax 



+ 0{e) 



where we have used the facts that p" does not depend on v and that L — and hence L ^ 
too — only applies to functions of v. Then, using this relation in (1^ . we get 



P^ 



with , 1 



At 

o". I 1 ' 



1 

Ax 



P 



i+i 



P'l 



A ^i-^' A 

Ax 2 /\x 

which is the usual 3-points stencil explicit scheme for the diffusion 



Pi Pi-i 



-(^A,iPi 



Si, 



equation (j4|). 

Of course, this property is true only if At can be chosen independently of e, or in other 
words, if the scheme is uniformly stable with respect to e. In [23], we have proved that 
this scheme is indeed uniformly stable under some CFL condition, in the simpler case of the 
telegraph equation (in which we have only two discrete velocities v = ±1, and s = 1, a = 1, 
aA = S = 0). In the following section, we extend this result to the general equation ([1]) for 
scheme (El). 
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3 Uniform stability 



We give our main result of stability for scheme ([6]). Without lost of generality, we assume 
S = (source free case). 

Theorem 3.1. If At satisfies the following CFL condition 

At<\ (aAx^ + 2eAx) , (9) 

with a = 2smO'm, then the sequences and g"' defined by scheme fM) satisfy the energy 
estimate 



7° 



Ax 



for every n, and hence the scheme ^ is stable. 

Note that CFL condition ([9]) can be viewed as an average of a diffusive CFL condition 
At < aAx"^ (needed for the diffusion scheme (IE])) and of a convection CFL At < sAx. It 
shows that the scheme is stable uniformly in e, that is to say a diffusive CFL condition 
At < CAx"^ is sufficient for stability for small £, while a convection CFL is sufficient for 
e = 0{l). 

Moreover, while our CFL condition ([9]) is valid for every e, we point out that it is not 
optimal for each e. This might be the price to be paid for getting a uniform condition. For 
instance, in the simpler case of constant cross sections (cr = (Tm) and kernel (s = Sm = |), the 
diffusion coefficient of the diffusion equation dH]) is k = 3^, and the optimal CFL condition 
for scheme (IE]) is At < However, CFL condition ^ for our AP scheme reads (for 5 = 0) 
At < 1^^, which is 0.2 as small as the optimal CFL. 

Remark 3.1. If we have a source term S 0, then the same CFL condition naturally gives 
a linear growth of the energy. Since this is standard and does not lead to any additional 
difficulty, we do not consider this case here. 

Remark 3.2. In the case of a time explicit discretization of the absorption term (that is 
to say when —crA,iP^~^^ is replaced by —(JA,iPi in (jSj), we have the same result with a CFL 
condition which is a bit more complicated. Namely, the scheme is stable if 



At < min — , — At5 , (10) 

where At^ is the maximum time step allowed by CFL condition (Q for the scheme with 
implicit (or zero) absorption term. 

This theorem is proved in the following sections: section 13.11 contains compact notations 
and useful lemma, and section 13.21 contains the derivation of our energy estimate. 
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3.1 Notations and useful lemma 

We give some useful notations for norms and inner products that are used in our analysis. 
For every grid function /i = (/ii)igz we define: 

ll/zf = 5^/i^Ax. (11) 
For every velocity dependent grid function v G [—1, 1] ^— 0(f) = (0j+i (f ))igz, we define: 

lll'^lll = EKi]^^- (12) 

If (/) and -0 are two velocity dependent grid functions, we define their inner product: 

(0,^) = 5^[0,^.i^,+ i] Ax. (13) 

Now we give some notations for the finite difference operators that are used in scheme ([6]). 
For every grid function = (0j4.i)jgZ5 we define the following one-sided operators: 

0i+l - 0i-l , 0i+3 - 0j . 1 

1 = ^\ ^ and 1 = ^\ (14) 

We also define the following centered operators: 

Finally, for every grid function \i = we define the following centered operator: 

^ = (16) 

Lemma 3.1 (Centered form of the upwind operator). For every grid function (p = {4>i^i.)iei! 
we have: 

{v+D- + V-D+) ct>,^. = vD^cP,^ - ^\v\D-D^<p,_,^. 

Proof. This result is easily obtained by using the relations = ^^^^ and the identity 
D- + = 2D^. □ 

Lemma 3.2. For every grid function (j) = (0j+i)igZ; we have: 

2 . 4 



J:{d^<P.^) Ax <^ 5^0^^. Ax. 
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Proof. Expand the left-hand side, use the inequahty (a + 6)^ < 2(a^ + b"^), and then use the 

2 ' ^ 2- 



change of index i + ^ ^ i + □ 



Lemma 3.3 (Estimate for the adjoint upwind operator). For every positive real number a 
and for every velocity dependent grid functions (p and if), we have: 

Proof. From Young's inequahty, we obtain for any positive real number a: 

\{{v+ D+ ^ V- D-) , 0)1 < a|||0||p + ^111 {v+ D+ ^ V- D-) . (17) 
Now, the second term of the right-hand side of this inequality can be written: 

^11 1 [v^D^ + V-D-) n? = ^Y.\{f AD-A^^^r dv + f v\D^^,^.f dv) Ax. 

Then, with a simple changing of index, D" can be replaced by in the first integral, which 
gives: 



^ ■■■\v\D-^ 



Finally, using this inequality in f|T7j) gives the result. □ 

Lemma 3.4 (Discrete integration by parts). For every grid functions (p = (0j+i)igZ; ^ = 
(V'i+i)iez, and fj, = {fii)i^z, we have: 

Proof. These results are simply obtained by using obvious changing of indexes and the 
definitions of the finite difference operators given in (ITH - fTBl) . □ 



Lemma 3.5. If g e L^([-l,l]) then 



[vgf < \ [\v\g'] 



8 




3.2 Energy estimates 

Since the time discretization of absorption term —a^p is implicit, it plays no role in the 
energy estimate. Consequently, to simplify the proof, we consider the case where there is no 
absorption (see remark [3731 at the end of this section). We proceed in five short steps. 
Step 1. 

Here we derive a first energy relation. With the finite difference operators defined in fll4p - 
(fTB]l . the scheme can be written in the following compact form: 



[vgr'] = 



At 



£2 £2 



(18a) 
(18b) 



We define the energy of system as dx + f [g'^] dx. It is clear that the scheme 
can be proved to be stable if the energy at time n + 1 can be controlled by the energy at 
time n. Consequently, we first multiply fll8al) by p""^^, then we take the sum over i G Z, and 
finally, we use the standard equality a(a — h) = \{a? — b'^ + \a — 6p) to get: 



2At 



|^n+i||2 _ ||^n||2 ^ ||^„+i _ ^„||2) + J] p^^D^ [vg^'] Ax = 0. (19a) 



Second, we multiply (118b|) by g^^i, we take the velocity average, we sum over i G Z, and we 



get: 



2At 



f+i II |2 - II |^"|| |2 + II |^"+i - |2) + _ (^"+1 , (/ _ [.]) (y-^D- + g-) 



6'pl.Ax. 



(19b) 



Now we use relation ([7]): since 



for every a simple expansion of the inner product 
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of the left-hand side of (119bp shows that it can be reduced to: 

{g'^+\ {I - []) {v+ D- + V- D+) g-) 



i&Z 

{g-^' , {v^D- + V-D+) g-) . 



Ax 



(20) 



Moreover, we can use relation and the assumptions of the cross section and the kernel to 
estimate the inner product of the right-hand side of fll9b[) as follows: 



Ax 



— '^^m^m II IQ 



n+l\\ |2 



(21) 



Consequently, we add up fll9ap and times (llQbp . then we use (l20l) . (12T!) . and the discrete 
integration by parts of lemma 13.41 to get our preliminary energy estimate: 



2At 



IP 



n+l ||2 



ip"f + iip"^^ - Pin + K"^'] ^ 



a; 



+ 



n+l|| |2 



2At 

< -5"|||fi'' 



n+l|| |2 



,"11 |2 



1^"+^ - <?"|| n + e {g-^' , + D^) g^) (22) 

5^ [vD'gr'] Pl^x. 



Step 2. 

In this step, we show how the p^^^ — p" term can be eliminated in (l22l) . First, it is use- 
ful to write p" in the right-hand side of (!22l) as (p" — p"^^) + Pi^"*"^: indeed, the terms 
Y.i&Pl^^D^ [vg^^^] and ^.^^ [fl^°^r^^] P?^^ in the left and right-hand sides cancel out 
and we obtain: 

|p«+l||2_||pn||2+||pn+l_pn||2) 



2At 



e2 



+ ^ (II \9^^'\\ ? - II l^?"ll r + II - 9l\?) + e (^?"+^ , iv^D- + t;-Z^+) (23) 

= -^iii^-^ir + E [^D'ar'] (pr - pr^)Ax. 

Now, we use the following Young inequality: 

E [^D'ar'] (Pr - Pr^) Aa: < allp-^^ - P"f + ^ E [^D'ar'] " Ax. (24) 



10 



Then the p"""*"^ — terms cancel out in (123|) if a = and we get 



1 



IP 



n+l ||2 



2At ^ 



iip"ir) 



II |2 



2At 



At 



Ax. 



(25) 



Step 3. 

Here, we work on the inner product of (!25|) to show that the g"''^^ — g"' terms can also be 
eliminated. First, we insert g'^'^^ in this inner product to get: 

(^"+\ {v+D- +V-D+) c/") 

= {g''+^ , {v+D~ + V-D+) ^"+^) + (c/"+^ , {v+D- + (c/" - ^"+^)) (26) 

= A + S, 

and we rewrite terms A and B as follows. For A, we use the centered form of the upwind 
operator (lemma [3.11) and the discrete integration by parts of lemma [33] to get: 



^{D^g-^\\v\D^g-^') 
Ax X - 

~9~ 2-^ 



Ax. 



For B, we also use the discrete integration by parts of lemma [231 to get: 

B = -{{v^D+ + V-D-) g''+^ , g'' - 6/"+^) . 
Then we apply the inequality of lemma 13.31 to B to get 



(27) 



(2^ 



\B\<aM'^'-gm' + —\\\\v\D^g-^'\\\'. 

4a 



(29) 



Therefore, using (125!) . (126!) . (127]) and (!29|) . we see that the ^r""^^ — ^r" terms cancel out in ([251) 



if a = and we get 



2At 



2At 



+4e h(«^c1 



n+l|||2 l||^^^|||2^ 



At, 



Ax - =^||||t;|DV^il^ 



(30) 



<-^iii/'^lir + 



At 



Ax. 
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Step 4- 

Now, we show how all the D^g"'^^ and the D^g'^^^ terms can be controlled by || |- 
First, note that the term ^|| | j-ylD+^i^+^ll p of the left-hand side of ([30]) can be estimated as 
follows: 



< 



Ax 
Ax, 



(31) 



since |f| < 1. Moreover, using lemma [375] and a change of indices shows that the last term 
of the right-hand side of (!30l) satisfies 



At 

T 



Ax. 



(32) 



Finally, we use these two estimates in (130!) to obtain: 



2At 



ip"^'ir-iip"f) + 



n+l|||2 |||„«|||2 



< -(^\\\9 



n+1 II |2 



2At 
/3At Ax 



~9' I 2-^ 



v\[D^g:^! 



Ax. 



(33) 



Now, taking the positive part of the factor (^ — e^) of the right-hand side of we 



have the estimate 
/3Af Ax\ sr^ 



Ax < 



< 



where we have used | f | < 1 and the estimate of lemma I3.2[ 
Step 5. 

Finally, estimates (133!) and (!34|) show that 



3At 


Ax\ 


+ r 

E 


4 










3At 




+ 4 


4 




Ax2' 



^n+l|| 1 2 



Ax 



(34) 



1 



IP' 



n+l ||2 II „n||2 



1^ 



n+l|||2 |||„n|||2 



2At ^'"^ " '"^ " ^ 2At 
This means that we have the final energy estimate 



/3At Ax 



Ax^ 



^ 111^ 



|p«+l||2 + £2|||^n+l|||2< ||^n||2^^2|||^n|||2 



n+l II |2 
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if At is such that 

'3 At Ax 



-e — ^-7 < a. 



4 2 J Ax 

Since 5" > 0, an equivalent condition is — < cr, which gives the sufficient 

condition 

Ax^a 2 , 

At < h -sAx, 

- 3 3 ' 

which proves the theorem. 

Remark 3.3. As explained at the beginning of this section, when the absorption term is 
non zero and is discretized imphcitly, its contribution —J2iez'^AiiP7~^^y^^ to energy 
estimate is non-positive and plays no role in the previous analysis. However, if we use instead 
an explicit discretization, then our analysis has to be modified. The difference now is that 
there is the additional term — ^jg^ o"A,jp"'''^p"Ax in the right-hand side of fl22l) . The idea 
of step 2 can be applied to this term (replace by — p""*"^) + p""*"^ and use a Young 
inequahty) so that the p"+^ — p"' terms cancel out, but now with a = 2At(i+o-A m) ' '^^^ other 
steps are the same, except that some coefficients are different. In order to shorten the paper, 
these details are left to the reader. 



4 Error estimates 

In this section, we simplify the presentation by taking constant total cross section a and 
kernel s = |, no absorption {a a = 0), and no source term (S = 0). These assumptions are 
not restrictive at all, and our analysis could be directly applied in the general case. 

Let T > be some finite time. For this study, we assume that the exact solution (p, g) 
of ([SD has the following regularity: 

peC\[0,TlH\R))nC\[0,T],H%R)) 

g E C\[0, T], IIH\R))) n C°([0, T],L\[-1., 1], ^ ' 

Note that this assumption implies that g is uniformly bounded with respect to e, which is 
quite strong. In particular, this excludes the case of initial or transition layers. Indeed, if, 
for instance, / is not isotropic at t = 0, then p{t = 0) 7^ f(t = 0), and hence g{t = 0) = 
-(/ — p)\t=o = 0{e~^) cannot be uniformly bounded with respect to e. 
With this assumption, we can obtain the following result. 

Theorem 4.1. If At satisfies the following condition 

AxV 2 , 

At < —— + -sAx, 36 
3 
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then for every time T > 0, there exists a constant C independent of At, Ax, and e, such 
that the numerical solution obtained by scheme ^ satisfies the following error estimate 



max 

n,nAt<T 



|Ax + £^ \ \g{tn, 



Ax 



< C((l + e2)At + Ax^ + eAx) 



This theorem is proved in the following sections: in section 14.11 we first derive the trun- 
cation error of the scheme, then in section IT2l we apply the same analysis as for the stability 
result to prove the theorem. 



4.1 Truncation error 



Let a" and ^6" be the truncation errors of scheme ([6]), that is to say the reminders obtained 
by inserting the exact solution of ([5]) in relations ([6]): 



p{tn+l,Xi) - p{tn,Xi) 

At 



g{tn+l,Xi^l) - g{tn+l,X^_i) 
J i — 

Ax 



(37a) 



At 



(37b) 



+ ^^(^~ [•]) {^^{9{tn,x^^^_,v) - g{tn,Xi_i,v) +v {g{tn,Xi^,v) - g{tn,Xi^,v))^ 



(37c) 



-^g{tn+i,Xi^i,v) - -^v 



1 p(t„,Xi+i) - p{tn,Xi) 1 



—6". 

2 « 



e^' ' - ' 2 ■ Ax 
We can prove the following estimate of these truncation errors: 
Lemma 4.1. There exists a constant C independent of At, Ax and e, such that 

||a"|| + II I < C ((1 + e^)At + Ax^ + sAx) 

for every n 

This lemma is proved by using standard simple techniques (Taylor-Lagrange formula of 
different orders). But to simplify the paper, the proof — which is a bit long — is given in 
appendix |Al 

Now, let p" and g^_^i be the convergence errors, that is to say the sequences defined by 



Pi = p{tn, Xi) - Pi and 1 (v) = g{tn, x^^ i , v) - g^^i (v) 
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Then these sequences satisfy the "perturbed" scheme, written in the following compact form: 



At 

~n+l ~n 

9i+i - 9i+i 1 



(38a) 



At 



+ ^(/ - [.]) [v^D- + v-D^) = -^g^^l - j,v5'~p-^^. + ^^h^, (38b) 



2 £^ 



with the homogeneous initial data = and g.^i = for every i. 



4.2 Analysis of the convergence error 

In this section, we apply the same analysis as for the stability result to prove that and g 
can be controlled by the truncation errors. 
Step 1. 

We multiply f l38al) by p""*"^ and take the sum over i to get 
1 



2At 



~n+l||2 _ II ~n||2 ^ ||~n+l _ ~n||2^) ^_ ^ ^^~n+lj _ ^ p^+^a^Ax. (39a) 



Second, we multiply (138bp by g""^} , we take the velocity average, we sum over i, which yields 



1 - '^"+1 II |2 _ II l^'^ll |2 + II |^"+1 _ ^"11 |2) + i <^n+l ^ _ f^]) ^ 
1 



2At 



--^111^ 



r,n+l\\ |2 



Finally, we add up fl39ap and times f l39bp , we use (120|) and lemma 13.41 to get 
1 



2At 



|^n+l||2 _ II ~n||2 ^ || _ ~n||2) + ^^"+1^0 ^^~n+lj ^ 



(39b) 



+ 



|~n+l|||2 |||;;n|||2 , \\\7,n+l 7,n\\\2 



2At 



Here, we can copy the steps 2 to 4 of the stability analysis (see section |3]) for relation ( 140|) . 
Therefore, skipping the details, we just give the resulting energy estimate: 



1 /n^„+l||2_||~n||2^)^ ^ 



2At 



|^n+l|||2 lll~ri|||2^ 



1 2 



+ 



2At 



Ax 



(41) 



+ j2P7<^^ + {r^\b-). 
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step 3. 

We estimate the scalar products of the right-hand side of fHTl) by using two different Young 
inequahties: 



(42) 



^pXAx<^(||p"+^|P+||a"|r) 

{r^\h-)<^-\w^'\\\' + ^\rr. 

Moreover, by taking the positive part of the factor — s^) in (HT]) . we have the estimate 



Ax < 



4 2 Ax2'"^ 



(43) 



Now, we can use the inequahties ( H2l) and ( H3l) to obtain the energy estimate 



< (1 + At)||p'^||2 + £^|||^"|||2+ -- + 



4 /3At Ax 



2' A.H4 l^A^lll^"^^'"^ 



(44) 



+ At||a"f + — I 

0" 



un\\ |2 



Note that the factor of || P in the right-hand side of (jHj) is non-positive if (^ — 

< |. Since u > 0, an equivalent condition is (^ — < |, which gives the sufficient 
condition 

AxV 2 ^ 

At < h -eAx, 

- 6 3 ' 

which is a bit more restrictive than the CFL condition ([9]). In that case, the energy esti- 
mate (jH]) can be simplified in 



\\~P 



^«+l||2 + ^2|||~n+l|||2 



<{l + At){\\~p-r + e'\\\r\\?)+^tC„{\\a 

where = 1 + -. 

Now, by using a simple recursion, we obtain 

iip"ir+^'iii^"iir 
<{i+/^tr{\\~pr+e'\m?) 



mi 2 



un\\ |2 



k\\2 



un-k\\ |2 



) (1 + Ai) 



fe-i 



fc=i 
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If n is such that nAt < T, we can use the classical inequality (1 + At)" < e"^* < e^. 
Moreover, we can use lemma W7\] and the fact that fP = = to get 

||p"f + £2|||^«|||2 

< e^C^TC ((1 + e^)At + Ax^ + sAxY , 

and hence 

llp"ll 

< C {{1 + e^)At + Ax^ + sAx) , 

where C = \fCJTC is independent of At, Ax, and e. This concludes the proof of the 
theorem. 

5 Conclusion 

In this paper, we have proposed a very simple stability proof for the recent AP scheme 
of • An explicit CFL condition that can be used in a computational code has been found: 
it insures that the scheme is stable and accurate, independently of e. This condition gives 
a classical parabolic CFL condition At = O(Ax^) when e goes to 0. Our proof uses very 
basic and simple arguments (energy estimates and Young inequalities), and is valid for one- 
dimensional linear equations with non-constant coefficients, continuous velocity variable, in 
the whole space. Our technique applies to general linear collision operators like operators of 
neutron transport or linear radiative transfer. By the same technique, we have also proved 
uniform error estimates. We mention that, in a work in preparation [23], we are able to 
apply our method to a simple non-linear problem coming from radiative transfer. 

In the future, the analysis of this scheme for initial boundary-value problems will be 
investigated. It would also be important to extend the scheme and its analysis to 2 or 3 
dimensional problems. 
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A Estimate of the truncation errors: proof of lemma [47T 



First, we remind (with no proof) that standard Taylor expansions give the following estimates 
for some time finite differences. 

Lemma A.l. (%) If ^ e C\[0,T]), then 

\ijitn+i) - i^itn)] < Atmax|^'|. 



Ifije C^[0,T]), then 



i^'itn)] < At max I V'" I . 



At V /I - jg^^j 

Then we also have the following estimates for some space finite differences. 
Lemma A. 2. (i) If (p e i/^(R) and Ax < 1, then 

5^0(x,)2Ax< 211011^.. 
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(a) If(t)e H'^{m), then 



E 



MXj+i) - (t){Xi) 

' Ax 



L2- 



(in) If(l)e H^{R), then 

E 



MXj+i) - (j){Xi) 

Ax 



-<i>'{x,^0\'Ax<^\\r\\h. 



Proof. For (i) , we write the difference of the continuous and discrete norms of as 
/ (l)(x)^ dx -'^(l){xifAx 



{ f2<^'{y)(t>{y)dy]dx. 



Then, using a simple Young inequahty, we get: 



J2 <l>{xi?Ax < / dx + J2 I (I {^\yf + <t>{yf) dy ) dx 

< / <j>{xr dx+j2 i {^\yf + dy) dx 

= ||0||i. + Ax(||0'||i. + ||<^||i.) 

<{\+Axmf^.. 

For (ii) , we use the Taylor-Lagrange formula up to first order to get 



<\>{xi^^) - <\>{xi) . 1 



Ax 



(f>'{Xi) = 



Ax 



{x — Xi+i)(f)"{x) dx. 



Then a simple Cauchy-Schwarz inequality gives the following: 



(f>{Xi+i) - (f){Xi) . 



Ax 



At 

rixfdx. 



Finally, we get the result by multiplying by Ax and taking the sum over i e Z. 
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For (iii), we use the Taylor-Lagrange formula up to second order to get the following two 
relations: 



{Xi) - (j){x 



""'+2 {Xj-XiY_ 

2Ax 



4>"'{x) dx. 



Then, we take the difference of these relations to obtain: 



(Xi+i) - (j){Xi) 





Ax 



»+2 



- x.^.r .... ^ r-i - ^0 



2Ax 



2Ax 



(t)"'{x) dx. 



Now, using a Young then a Cauchy-Schwarz inequalities gives 



(Xj+i) - (j){Xi) 



Ax 



- 9 \Xi+\ 



< 



Ax^ 
320 



3 r^i+i 



<P"'{xYdx. 



Again, the final result is obtained by multiplying by Ax and taking the sum over i Eli. □ 



Now we study the truncation error a" which is defined by ( 137ap . Since (p, g) is the exact 
solution, we have: 



p{tn+l,Xi) - p{tn,Xi) 

At 



dtp{tn,Xi) + 



g{tn+l,Xi^l) - g{tn+uXi_i) 



Ax 



Then we estimate the norm of a" as follows: we use a Young inequality, lemma IA.lt (iii) of 
lemma IA.2t and lemma 13.51 and we obtain 



^KpAx < 2 I ^ At2max|atip(t,Xi)|^ Ax + ^ 



Ax^ 
320" 



\dxxxg{tn+l, ■,V) 



lL2 



By using (i) of lemma IA.2t the first term of the right-hand side of the previous estimate can 
be controlled by a norm of p (if Ax < 1) and we get 



5]]K|'Ax < 2 (2At^\\pfx + 

ic7 V 



1 Ax^ 



4 320 



where X = C\[0,T], H\R)) and y = C^{[0,T], L\[-l,l], H%R))). Consequently, it is 
clear that ||a"|| < Ci{At + Ax^), where C*i depends only on the norms of p and g. 
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Now, we study the truncation error 6" which is defined by f l37cp . Again, since {p,g) is 
the exact solution, we have: 

&r = £ I dtgitn,Xi,v) 

+ e{I-[.]) I v+ \ ^ — ^ d^g{tn,Xi_i,v) 

_ g{t„,Xi^,v)-g{tn,Xi^i,v) 
+v I — d,g{t^,x^^,v 

+ (gitn+l,Xi^l,v) - g{tn,Xi^,v) 

[ d,p{tn,x,^ 

The sequel is similar, though a but longer, to what we did for a". By using again lemmas [A. II 
and IA.2t and some standard inequalities (Young and Jensen), the reader can easily find that 
the following estimate is satisfied: 

Y: mn Ax < 4 (.WH.r^ + £^2^Ax1,||^ + aWH^H^ + ^IIpII^) , 



where X = C^{[0,T], L'^{[-1,1], H\R))) and Y = C"{[0,T], H^{R)). Therefore, we have 
II 1^*^11 1 ^ C2 {{1 + e^)At + eAx + Ax"^), where C2 depends only on the norms of p and g. 
Finally, the estimates of ||a"|| and |||b"||| allow us to complete the proof. 
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